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Abstract.  The  purpose  of  this  paper  is  to  illustrate  that  direct  simulation  Monte  Carlo  methods  can  often  he  considered  as 
rigorous  mathematical  tools  for  solving  nonlinear  kinetic  equations  numerically.  First  a  convergence  result  for  Bird’s  DSMC 
method  is  recalled.  Then  some  sketch  of  the  history  of  stochastic  models  related  to  rarehed  gas  dynamics  is  given.  The  model 
introduced  by  Leontovich  in  1935  provides  the  basis  for  a  rigorous  derivation  of  the  Boltzmann  equation  from  a  stochastic 
particle  system.  The  last  part  of  the  paper  is  concerned  with  some  recent  directions  of  study  in  the  field  of  Monte  Carlo 
methods  for  nonlinear  kinetic  equations.  Models  with  general  particle  interactions  and  the  corresponding  limiting  equations 
are  discussed  in  some  detail.  In  particular,  these  models  cover  rarehed  granular  gases  (inelastic  Boltzmann  equation)  and  ideal 
quantum  gases  (Uehling-Uhlenbeck-Boltzmann  equation).  Problems  related  to  the  order  of  convergence,  to  the  approximation 
of  the  steady  state  solution,  and  to  variance  reduction  are  briehy  mentioned. 


DSMC  AND  THE  BOLTZMANN  EQUATION 

Direct  simulation  Monte  Carlo  (DSMC)  has  been  the  most  widely  used  numerical  algorithm  in  kinetic  theory.  We  refer 
to  G.A.  Bird’s  monograph  [1,  Sections  9.4,  11.1]  concerning  remarks  on  the  historical  development.  The  history  of  the 
subject  is  also  well  reflected  in  the  proceedings  of  the  bi-annual  conferences  on  “Rarefied  Gas  Dynamics”  ranging  from 
[2]  to  the  present  volume. 

The  method  is  based  on  systems  of  particles 

t>0,  n>l,  (1) 

imitating  the  behaviour  of  gas  molecules  in  a  probabilistic  way.  It  includes  several  numerically  motivated  approxima¬ 
tions.  Independent  motion  (free  flow)  of  the  particles  and  their  pairwise  interactions  (collisions)  are  separated  using  a 
splitting  procedure  with  a  time  increment  At .  During  the  free  flow  step,  particles  move  according  to  their  velocities, 

/t  +  At 

Vi{s)ds,  i  =  l,...,n, 

and  do  not  collide.  At  this  step  boundary  conditions  are  taken  into  account.  During  the  collision  step,  particles  do  not 
change  their  positions.  At  this  step  some  partition 

of  the  spatial  domain  into  a  finite  number  Ic  of  disjoint  cells  is  introduced.  In  each  cell  a  certain  amount  of  binary 
collisions  between  particles  is  performed, 

Vi,Vj  ^  V*,V*. 

The  probabilistic  rules  of  these  transformations  depend  on  the  interaction  potential  between  gas  molecules. 

The  interest  in  studying  the  connection  between  stochastic  simulation  procedures  in  rarefied  gas  dynamics  and  the 
Boltzmann  equation  was  stimulated  by  K.  Nanbu’s  paper  [3]  (cf.  the  survey  papers  [4,  5]).  Convergence  for  the  Nanbu 
scheme  and  its  modifications  was  studied  in  [6]  (spatially  homogeneous  case)  and  [7]  (spatially  inhomogeneous  case). 
Convergence  for  Bird’s  method  was  proved  in  [8].  It  was  established  that  the  empirical  measures  related  to  the  particle 
system  (1), 

1  " 

p,^'^\t,dx,dv)  =  -'^5^^(t),vi(t){dx,dv) , 
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converge  (as  n  oo)  to  a  deterministic  measure  with  some  density  /,  which  is  defined  as  follows. 

Consider  the  time  discretization  tk  =  k  At ,  k  =  0,1, . . . ,  and  auxiliary  functions 

,  te[tk,tk+i],  {x,v)eDxn^. 

These  functions  are  determined  by  two  systems  of  equations  coupled  via  their  initial  conditions.  The  first  system, 
corresponding  to  the  free  flow  simulation  steps,  has  the  form 

^  X,  v)  +  {v,  x,v)  =  0,  (2) 

with  appropriate  boundary  conditions.  Here  V  denotes  the  vector  of  partial  derivatives  and  ( ., . )  is  the  scalar  product 
in  the  Euclidean  space  TZ^ .  The  initial  conditions  are 

f^^’'"'>{tk,x,v)  =  f^^’’^~'^'>{tk,x,v),  k=l,2,...,  f^^’°'>{0,x,v)  =  fo{x,v). 


where  /o  is  some  density  function  related  to  the  initial  state  of  the  particle  system  (1).  The  second  system  of  equations, 
corresponding  to  the  collision  simulation  steps,  has  the  form 


d 

dt 


f^^’'^\t,x,v)  = 


where 


(3) 

f^^’’^'>{t,x,v)  f^^’’^'>{t,y,w)^dedwdy. 


h{x,y) 


1 

W\ 


Ic 

J2^Di{x)xDi{y), 

1=1 


(4) 


I .  I  denotes  the  volume,  xa  is  the  indicator  function  of  a  set  A  and  is  the  unit  sphere  in  TZ^ .  The  collision  kernel  B 
is  determined  by  the  interaction  potential  between  gas  molecules.  The  initial  conditions  are 

(ffc,  v)  =  (ffc+i,  x,v),  k  =  0,l,.... 

The  limiting  density  /  is  defined  as 

f{tk,x,v)  =  f^^'’"~^\tk,x,v),  k=l,2,...,  f{0,x,v)  =  fo{x,v). 

HISTORY  OF  STOCHASTIC  MODELS  IN  RGD 


Statistical  description 

The  basis  for  the  statistical  theory  of  gases  was  provided  in  the  second  half  of  the  19th  century.  James  Clerk  Maxwell 
(1831-1879)  found  the  distribution  function  of  the  gas  molecule  velocities  in  thermal  equilibrium. 


feq{v)  = 


1 


[27rT]3/: 


:  exp  - 


2T 


veTZ^ 


(5) 


The  positive  number  T  corresponds  to  the  temperature  of  the  gas.  Ludwig  Boltzmann  (1844-1906)  studied  the  problem 
if  a  gas  starting  from  any  initial  state  reaches  the  Maxwellian  distribution  (5).  In  [9]  he  established  the  equation 


d 

dt 


f{t,x,v)+(v,V,^)f{t,x,v)  =  J  j  ^deB{v,w,e)\j{t,x,v*)  f{t,x,w*)-f{t,x,v)  f{t,x,w)\^  (6) 


governing  the  time  evolution  of  the  distribution  function 


t>0,  X  e  D  cTZ^ ,  V  £TZ^  . 
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This  function  represents  the  relative  amount  of  gas  particles  at  time  t  having  a  position  close  to  x  and  a  velocity  close 
to  V .  The  collision  transformation  is 

V*  =  V  +  e{e,w  —  v) ,  w*  =  w  +  e(e,v  —  w) .  (7) 

For  the  hard  sphere  model  (billiard  balls)  the  collision  kernel  takes  the  form 

B[v,w,e)  =  c\{v  —  w,e)\ ,  for  some  c>0. 


Leontovich  model 

Consider  a  stochastic  particle  system 

zW(f)=  14(f)),  f>0,  (8) 

determined  by  an  infinitesimal  generator  of  the  form 


n 

A^'^\<^){z)  =  Y,{vi,V,,m{z)  + 

—  E  / 

2  77  1  ^ 

t)) 

-$(2)j  q^'^'>{xi,Vi,Xj,Vj,e)de 

i  =  l 

where  $  is  an  appropriate  test  function. 

2  =  {xi,Vl-,.. 

•  5  ^n)  5 

i  = 

is  an  intensity  function  and 

f  [xk.Vk) 

5 

if 

k  i,j, 

[J{z,i,j,e)]k  = 

{xi,Vi  +  e{e, 

Vj  -Vi))  , 

if 

k  =  i, 

y  {xj,Vj-fe{e, 

.Vi-Vj)), 

if 

k  =  j, 

is  the  jump  transformation.  Here  we  assume  D  =  TZ^  so  that  no  boundary  conditions  are  involved.  We  recall  the 
path  wise  behaviour  of  the  system  (8),  which  is  a  piecewise-deterministic  process  in  the  sense  of  [  1 0,  Section  2] .  Starting 
at  2  the  process  moves  according  to  the  free  flow,  i.e. 

^(f)  —  T f  r’l ,  , . . . ,  T f  ,  1^^)  , 


until  a  random  jump  time  ri  is  reached.  The  probability  distribution  of  this  time  is  determined  by 

ft 


where 


Prob(ri  >  f)  =  exp  y  ds^  ,  t>0, 


1  <  1 j  <  n 


At  the  random  time  ti  the  process  jumps  into  a  state  zi ,  which  is  obtained  from  the  state  z'  =  . .  .-jx'^jv!^)  of  the 

process  just  before  the  jump  by  a  two-particle  interaction.  Namely,  two  indices  i,  j  and  a  direction  vector  e  are  chosen 
according  to  the  probability  density 


q('^\x'„v'„x'.,v'.,e) 

2nA(")(2') 

and  the  velocities  ,  n'  are  replaced  using  the  collision  transformation  (7). 


(9) 
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The  study  of  the  relationship  between  the  process  (8)  and  the  Boltzmann  equation  (6)  was  started  by  M.  A.  Leontovich 
(1903-1981)in  the  paper  [11]  in  1935. Let  z)  denote  the  n-particle  distributionfunction  of  the  process.  Applying 

properties  of  the  collision  transformation  (7)  and  some  symmetry  assumptions  on  ,  the  equation 

+  =  (10) 

i  =  l 

<n  ^ 

is  obtained  from  Kolmogorov’s  forward  equation.  Using  the  method  of  generating  functions,  Leontovich  first  studied 
the  cases  of  “monomolecular  processes”  (independent  particles)  and  of  “bimolecular  processes”  with  discrete  states 
(e.g.,  a  finite  number  of  velocities).  Under  some  assumptions  on  the  initial  state,  he  showed  that  the  expectations 
of  the  relative  numbers  of  particles  in  the  bimolecular  scheme  asymptotically  (as  n  oo)  solve  the  corresponding 
deterministic  equation.  In  the  case  of  the  full  Boltzmann  equation  the  stochastic  process  was  described  via  (10)  (even 
including  a  boundary  condition  of  specular  reflection).  Concerning  the  asymptotic  behaviour  of  the  process,  Leontovich 
noted  that  he  was  not  able  to  prove  a  limit  theorem  in  analogy  with  the  discrete  case,  though  he  strongly  believes  that 

(n) 

such  theorem  holds.  However,  he  pointed  out  the  following.  Let  p^.  '  denote  the  marginal  distributions  corresponding 
to  the  density  .  If 


lim  ®2,  *^2)  =  lim  »i,  «^i)  lim  »2,  *^2) 

n—^oo  n—^00  n—^00 


(11) 


and 


lim  q^'^\x,v,y,w,e)  =  5[x  —  y)  B[v,w,e) ,  (12) 

n— >-oo 

where  S  denotes  Dirac’s  delta-function,  then  the  function 

f{t,x,v)  =  lim  p^^\t,x,v) 

n—^oo 

solves  the  Boltzmann  equation.  For  example,  (12)  is  satisfied  for 

q^'^\x,v,y,w,e)  =  h^'^\x ,  y)  B {v ,  w ,  e) 


and 


if  ||®-2/||  < 
otherwise , 


where  c„  is  the  volume  of  a  ball  with  radius  0 .  In  this  case,  according  to  (9),  only  those  particles  can  collide 

which  are  closer  to  each  other  than  the  interaction  distance  . 

Independently,  the  stochastic  approach  to  the  Boltzmann  equation  was  developed  by  M.  Kac  (1914-1984)  in  the 
paper  [12]  in  1956  (see  also  the  books  [13]  and  [14]).  The  spatially  homogeneous  version  of  equation  (10)  is  called 
“master  equation”  referring  to  [15].  The  factorization  property  (11)  is  called  “chaos  property”  (indicating  asymptotic 
independence)  and  the  statement  of  the  basic  theorem  is  “propagation  of  chaos”.  The  spatially  inhomogeneous  case  was 
treated  by  C.  Cercignani  in  the  paper  [16]  in  1983.  He  considered  a  system  of  “soft  spheres”,  where  “molecules  collide 
at  distances  randomly  given  by  a  probability  distribution”,  and  proved  propagation  of  chaos  (modulo  a  uniqueness 
theorem).  The  limiting  equation  is  the  mollified  Boltzmann  equation  (cf.  [17,  Section  VIII. 3]) 


—  f{t,x,v)+{v,V^)f{t,x,v)  = 


(13) 


/  /  /  h{^,y)B{v,w,e)  f{t,x,v*)f{t,y,w*)-f{t,x,v)f{t,y,w)  dedwdy, 

Jd  Js^  L  J 


which  reduces  to  the  Boltzmann  equation  (6)  if  the  “mollifier”  /i  is  a  delta-function.  Further  references  concerning  the 
spatially  inhomogeneous  case  are  [18,  19,  20,  21,  22,  23,  24]. 
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Stochastic  numerics 


Stochastic  numerical  procedures  in  rarefied  gas  dynamics  go  back  to  the  papers  [25]  (homogeneous  gas  relaxation 
problem)  and  [26]  (shock  structure  problem)  by  G.A.  Bird  in  1963  and  1965,  respectively.  An  interesting  development 
of  the  DSMC  method  was  related  to  the  “time  counter”  -  the  mechanism  that  determines  the  number  of  collisions  during 
a  time  step.  Several  similar  modifications  were  introduced  under  different  names  in  [27]  (null-collision  technique), 
[28]  (majorant  frequency  scheme),  [29]  (no  time  counter  scheme).  The  basic  idea  of  “fictitious  jumps”  is  applicable  in 
general  situations  (cf.  [30,  Section  4.2]).  In  [3]  K.  Nanbu  derived  a  stochastic  particle  scheme  starting  directly  from  the 
Boltzmann  equation.  However,  the  original  method  suffered  from  certain  deficiencies  (quadratic  effort  in  the  number 
of  particles,  conservation  of  momentum  and  energy  only  on  average).  Later  this  approach  was  considerably  improved 
(cf.  [31,  6])  so  that  it  did  successfully  work  in  applications  like  the  reentry  problem  (cf.  [32,  33]).  Among  the  numerous 
Russian  sources  on  stochastic  algorithms  for  the  Boltzmann  equation  we  mention  papers  related  to  the  Leontovich- 
Kac -process  [34,  35,  36,  37],  papers  using  branching  processes  [38,  39]  and  papers  following  Skorokhod’s  approach 
via  stochastic  differential  equations  with  respect  to  Poisson  measures  [40,  41]. 


Comments 

Boltzmann’s  equation  (6)  was  derived  from  a  model  with  deterministic  interaction  between  particles.  The  origin  of 
stochasticity  (necessary  for  the  interpretation  of  the  solution)  remained  unclear.  Statistical  mechanics  solves  this  prob¬ 
lem  via  a  random  initial  state  of  the  system.  Leontovich’s  motivation  for  introducing  his  model  was  to  provide  a  stochas¬ 
tic  process  such  that  the  Boltzmann  equation  describes  the  limit  (for  the  particle  number  going  to  infinity)  of  rigorously 
defined  mathematical  expectations.  Moreover,  this  approach  allows  one  to  study  fluctuations.  He  did  not  address  the 
problem  of  the  relationship  between  the  stochastic  process  and  the  description  of  the  system  by  quantum  (or  classical) 
mechanics.  The  main  new  point  in  this  model  are  the  stochastic  interactions.  Billiard  balls  can  penetrate  each  other 
colliding  only  with  a  certain  probability.  This  is  a  conceptual  barrier,  which  Kac  was  not  able  to  overcome.  In  1959  he 
wrote  [14,  p.l31]:  “The  primary  disadvantage  of  the  master  equation  approach  ...  lies  in  the  difficulty  (if  not  impossi¬ 
bility!)  of  extending  it  to  the  nonspatially  uniform  case.”  In  1973  the  message  was  the  same  [42,  p.385]:  “The  master 
equation  approach  suffers  from  a  major  deficiency.  It  is  limited  to  the  spatially  homogeneous  case.  It  seems  impossible 
to  bring  in  streaming  terms  while  at  the  same  time  treating  collisions  as  random  events.”  Bird’s  DSMC  algorithm  can  be 
derived  from  the  Leontovich  model  by  introducing  certain  numerical  approximations  (splitting,  cells).  The  problem  of 
combining  free  flow  with  stochastic  collisions  is  avoided  (hidden)  by  the  splitting  procedure.  Nevertheless,  all  particles 
in  a  cell  collide  according  to  some  probabilistic  rules  and  not  according  to  their  actual  trajectories. 


SOME  RECENT  DIRECTIONS  OF  STUDY 
General  interactions 

Consider  particles  with  types  from  a  locally  compact  separable  metric  space  Z .  Any  event  in  the  system  consists  in 
the  interaction  of  at  most  R  particles  and  produces  as  a  result  a  finite  number  of  particles.  This  includes,  for  example, 
the  generation  of  new  particles  from  a  source,  the  extinction  or  transformation  of  single  particles,  the  collision  of  two 
particles,  as  well  as  more  complicated  chemical  reactions.  Formally,  we  introduce  the  state  space  of  the  particle  system 
as 


:  n>0,  Zi^Z,  i=l,...,nj, 

where  Sa  denotes  the  delta-measure  concentrated  at  a .  The  admissible  subsequent  states  of  zz  e  i?  are  denoted  by 

where  are  pairwise  distinct  indices  from  { 1 , . . . ,  n}  and  {  £  i? .  The  intensity  of  these  events  is  determined 

by  some  compactly  bounded  kernels  qo,qi, . .  .,qR  depending  on  a  corresponding  number  of  arguments .  The  asymptotic 
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behaviour  of  such  systems  has  been  studied  under  appropriate  assumptions  in  [43].  It  turned  out  that  they  approximate 
the  deterministic  measure-valued  solution  of  the  equation 


d 

dt 


=  t>o. 


(14) 


The  nonlinear  operator  is  defined  as 

R 


U-LL  [(<^,0  -  Pizi)  -  ...  -  )]  qr  [zi,.  ..,Zr,d^)  /i(d2i)  .  .  ./i(d2^)  , 

for  continuous  test  functions  (p  with  compact  support  and  measures  ii  on  Z ,  where  the  notation  {(p,  v)  =  p[z)  v[dz) 
is  used.  As  an  example  we  consider  the  spatially  homogeneous  Boltzmann  equation  with  dissipative  collisions.  Let 
Z  =  72.  ,d,>l,  and  denote 


V  {v,w,e,d)  =  — - l-s{v,w,d) - - - ,  w  {v,w,e,d)  =  — - s{v,w,d) - - - , 

where  v,w  ^  TV^ ,  e  £  (unit  sphere),  0  £  0  (measurable  space)  and  e  is  some  measurable  function.  Note  that 


y\?  +  \\w'\?  =  ||d|2  +  ||«;||2_ 


1  —  e[v,  w,  Oy 


i.e.  energy  is  dissipated  if  <  1 ,  conserved  if  =  1 ,  and  created  if  >  1 .  Assuming 

q2{v,w,A)  =  ^  J  J ^  ^  XAiSvi(v,.ui,e,e)  +  K'{v,w,e,e))B{v,w,e)deTY{de), 
where  A  is  a  measurable  subset  of  E  and  tt  is  some  probability  measure  on  0 ,  equation  (14)  takes  the  form 

—  /  /  /  \p[v' [v,w,e,6))  +  p[w' [v,w,e,6))  —  p[v)  —  p[w)\B[v,w,e)  deTT[d6)  ii[t,dv)  ii[t,dw) 

2  J-Rd  J-j^d  Jq  Jsd-1  L  J 

and  reduces  to  the  weak  form  of  the  spatially  homogeneous  Boltzmann  equation  in  the  case  e  =  1 . 

An  even  more  general  class  of  stochastic  models  is  obtained  when  the  intensity  of  an  individual  interaction  depends 
on  the  state  of  the  whole  system,  e.g.  via  the  particle  density.  Such  systems  approximate  certain  generalized  Boltzmann 
equations.  One  example  is  the  Uehling-Uhlenbeck-Boltzmann  equation  [44]  related  to  ideal  quantum  gases.  In  the 
spatially  homogeneous  case  this  equation  reads 

^  ^  deB{v,  w,  e)  ^[1+6 f{t,  ?;)]  [1+6 f{t,  w)]  f{t,  v*)f{t,  w*)  - 

[1+6 f{t,  ?;*)]  [1+6 f{t,  w*)]  f{t,  v)  f{t,  u;)j .  (15) 


It  reduces  to  the  Boltzmann  equation  in  the  case  0  =  0.  Bose-Einstein  statistics  and  Fermi-Dirac  statistics  are  equi¬ 
librium  distributions  of  equation  (15)  with  0  =  -|-1  and  0  =  —  1 ,  respectively.  A  DSMC  method  for  this  equation  was 
studied  in  [45],  where  more  details  and  further  references  can  be  found.  A  second  example  is  the  limiting  equation  for 
the  Consistent  Boltzmann  Algorithm  [46],  a  DSMC  modification  for  dense  gases.  This  equation  takes  the  form 


d 

—  f{t,x,v)+(v,Vyf{t,x,v)  = 


(16) 


[  dw  [  deB{v,w,e)  Y{g{t,x*))  f{t,x*  ,v*)  f{t,x*  ,w*) -Y{g{t,x))  f{t,x,v)  f{t,x,w)  . 

Jrs  Js^  l  J 


Here  g{t,x)  =  f[t,x,v)dv  is  the  density  and  x*  denotes  an  appropriate  displacement  during  collisions.  The 
positive  continuous  function  Y  is  going  to  infinity  as  the  density  reaches  its  close-packing  value.  Equation  (16),  which 
is  closely  related  to  the  Enskog  equation,  has  been  studied  in  [47,  48],  where  more  details  and  further  references  can  be 
found  (see  [49,  50]  concerning  DSMC  methods  related  to  the  Enskog  equation).  The  derivation  of  equations  (15)  and 
(16)  from  corresponding  particle  systems  was  rather  formal,  and  rigorous  results  would  be  of  considerable  interest. 
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Convergence  order 


The  DSMC  algorithm  mentioned  above  depends  on  three  main  approximation  parameters  -  the  number  of  particles 
n ,  the  splitting  time  step  At  and  the  cell  size  Ax .  By  taking  the  limit  n  ^  oo  one  obtains  the  system  of  equations  (2), 
(3),  which  depends  on  the  remaining  two  parameters.  In  the  limit  At  ^  0  this  system  transforms  into  equation  (13) 
with  h  defined  in  (4).  Finally,  the  Boltzmann  equation  (6)  is  obtained  in  the  limit  Ax  0 .  For  numerical  purposes 
the  order  of  convergence  with  respect  to  the  various  parameters  is  important.  The  order  with  respect  to  n  has  been 
established  in  the  context  of  general  Markov  processes  in  [20,  23].  Cell  size  dependence  was  studied  in  [51].  Second 
order  with  respect  to  time  step  was  predicted  in  [52] .  However,  in  [53, 54]  the  author  noticed  a  mistake  in  the  derivations 
of  [52]  and  gave  some  specific  example  illustrating  first  order  behaviour.  On  the  other  hand,  rather  thorough  numerical 
experiments  in  [55]  showed  second  order  behaviour  in  many  test  cases  (see  also  the  derivation  in  [56]).  This  problem 
was  further  investigated  in  [57].  See  also  the  contribution  by  M.A.  Gallis  and  co-authors  in  these  proceedings. 


Steady  state  problems 


In  many  applications  studying  the  equilibrium  behaviour  of  gas  flows  is  of  primary  interest.  To  this  end,  time  aver¬ 
aging  over  trajectories  of  the  corresponding  particle  system  is  used. 


1 

k 


k 


E 


i  =  l 


tj  =  t  +  j  At , 


where  (p  is  some  test  function  and  t  is  the  starting  time  for  averaging.  To  justify  this  procedure  (for  k  oo),  one  has 
to  study  the  connection  between  the  stationary  density  of  the  process  and  the  stationary  Boltzmann  equation.  From  the 
results  mentioned  above  one  can  obtain  information  about  the  limit  limt_yoo  lim„_yoo  x,  v)  while  here  one  is 

interested  in  the  limit  lim„_yoo  limt_yoo  [t^x^v) .  The  identity  of  both  quantities  is  not  at  all  obvious.  Some  result 
concerning  the  factorization  of  the  equilibrium  distribution  was  obtained  in  [58].  But  in  this  field  the  state  of  the  theory 
is  rather  unsatisfactory. 


Variance  reduction 

A  serious  problem  in  many  applications  of  DSMC  (flows  with  high  density  gradients,  or  low  Mach  number  flows)  are 
large  statistical  fluctuations.  For  the  purpose  of  variance  reduction,  a  modification  of  DSMC  called  SWPM  (stochastic 
weighted  particle  method)  was  proposed  in  [59].  In  SWPM  a  system  of  weighted  particles  is  used,  which  allows  one  to 
resolve  low  density  regions  with  a  moderate  number  of  simulation  particles  (cf.  [60]).  A  convergence  proof  for  SWPM 
has  recently  been  proposed  in  [61].  Further  references  related  to  weighted  particles  are  [62, 63].  An  “information  preser¬ 
vation  method”  (cf.  [64,  65])  has  been  developed  for  low  Mach  number  flows  occurring  in  micro-electro-mechanical 
systems  (MEMS).  Here  the  convergence  issue  is  open. 
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